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We study indented spherical colloids, interacting via depletion forces. These systems exhibit 
liquid-vapor phase transitions whose properties are determined by a combination of strong “lock- 
and-key” bonds and weaker non-specific interactions. As the propensity for lock-and-key binding 
increases, the critical point moves to significantly lower density, and the coexisting phases change 
their structure. In particular, the liquid phase is “porous”, exhibiting large percolating voids. The 
properties of this system depend strongly on the topological structure of an underlying bond network: 
we comment on the implications of this fact for the assembly of equilibrium states with controlled 
porous structures. 


In colloidal systems, a wide range of structures can 
be self-assembled from a simple palette of components 
and interactions [1]. For example, anisotropic attractive 
forces between particles can stabilise colloidal micelles or 
exotic crystals [2-7], which could form the basis of fu¬ 
ture materials and devices. Anisotropic attractions can 
be achieved experimentally via chemical patches [2, 6, 8- 
10], or through a combination of particle shape and de¬ 
pletion forces [4, 5, 11-13]. Depletion forces arise when 
colloids are dispersed in a solution of much smaller “de- 
pletant” particles, which induce an attractive interac¬ 
tion [14], whose properties can be tuned via the size and 
number density of the depletant. Here, we present Monte 
Carlo simulations of the assembly of indented (“lock-and- 
key”) colloidal particles [4, 15-18], in the presence of an 
ideal depletant [19]. We introduce an effective potential 
which captures quantitatively the effect of this depletant. 
This enables efficient simulation, and accurate character¬ 
isation of the liquid-vapor phase transitions that occur in 
this system. When the depletant particles are small in 
size compared to the colloids, the resulting liquid and 
vapor phases have unusual properties. In particular, the 
critical point for the phase transition occurs at a rather 
low density, and the liquid phases have a complex struc¬ 
ture that includes large voids. Indeed, we find that the 
liquid phase is “porous” on a length scale comparable 
with the collidal particles. Similar properties have also 
been found in model polymer liquids [20-22] and in the 
“empty liquids” that occur in patchy-particle models [23- 
25]. However, the situation for indented colloids is dif¬ 
ferent from both these cases: we trace these differences 
to the topology of the bond networks that appear as self- 
assembly takes place. 

An indented colloidal particle is modelled as a sphere 
of diameter a, from which we cut away its intersection 
with a second sphere of the same diameter [17, 18]. The 
distance between the sphere centres is taken as dc = 0.6cr, 
so the depth of the indentation (measured from the “lip”) 
is 0.2cr, comparable with experimentally realizable par¬ 
ticles [4]. The orientation of particle i is specified by a 
unit vector rii that points outwards through the centre of 
the indentation. The hard particle interactions between 
these particles are treated exactly in our model, but we 
parameterise the depletion interaction between the col¬ 



FIG. 1 . Illustration of the effective potential between indented 
colloids, (a) The interaction depends on the vector ri2 be¬ 
tween the colloids, and their orientations ni,2. We define two 
angles and 4 > by cos0r = niax(ni • ri2,—n2 • ri2) and 
coscf = rii • n2. Particles can bind in a lock-key configura¬ 
tion (small 0r) or a back-to-back configuration (large 0r), as 
indicated by the two shaded areas, (b) In the lock-key ori¬ 
entation, the interaction depends on the angle 4 > through the 
function W{(j>), the form of which is motivated by studies of 
the exact depletion potential (see Appendix A). 


loids through an effective potential, illustrated in Fig. 1. 
Specifically, when two colloids approach one another in 
the “back-to-back” configuration, they feel a square-well 
interaction potential of depth £bb and range cbb • When 
they approach in the “lock-key” configuration, the in¬ 
teraction range is clk, and the effective potential also 
depends on the angle (j) between the particle directors. 
For perfect lock-key binding, the depth of the potential 
well is £lk- We fix the unit of energy by setting /3 = 
l/{kBT) = 1. The interaction ranges are fixed through¬ 
out this work according to (clk —dc) = (t'bb —o') = O.ltr, 
comparable with the size of the depletant particles. (As 
for spherical particles [26], the behavior depends weakly 
on these interaction ranges, if the interaction strengths 
are adjusted so that the binding free energies or second 
virial coefficients for pair interactions are held constant, 
and the ranges are small compared to a.) To obtain 
the interaction strengths £lk and £bb as a function of 
depletant parameters, we used the geometric cluster al¬ 
gorithm (GCA) [27, 28] to simulate systems of two col- 
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FIG. 2. (a) Summary of phase behaviour, as a function of elKjEbb- In the shaded region, coexistence of liquid and vapor 
phases can occur, for suitable colloid densities p. The stars indicate the positions of critical points. The dotted lines show 
how the effective interaction strengths (slk, £bb) depend on the depletant size ratio q (the depletant volume fraction rja varies 
along these lines). (b,c) Representative configurations from the coexisting “liquid” and “vapor” phases at the indicated state 
point, which corresponds to Pb ~ 1.057;“. For the liquid state, the void space is colored blue, to emphasise the porous structure, 
(d) Binodal curves showing the densities of coexisting liquid and vapor phases as a function of depletant parameters q, pa and 
(total) colloid density p. Stars mark critical points. The dashed lines show the densities at which the void spaces percolate the 
system, as discussed in the main text. (e,f) Visualisations of the largest void spaces from representative configurations at the 
indicated state points, both of which correspond to t]b = 1.05r;s. 


loidal particles interacting with a depletant of penetrable 
spheres [19]. The size ratio between colloids and deple¬ 
tant particles is q (the depletant diameter is qa), and the 
depletant volume fraction is r]s. Full details of the effec¬ 
tive potential and the parameterisation of (elKj^bb) are 
given in Appendix A. 

We employed grand canonical Monte Carlo (GCMC) 
simulation to study the phase behaviour [29, 30] of the 
colloids, interacting through the effective potential, in 
cubic boxes of sizes between (12 (t)^ and (20cr)^. Smaller 
boxes sufficed for large g, where the size of bound clusters 
is typically small; bigger boxes were needed to accommo¬ 
date the larger clusters that form at small q. Our GCMC 
method uses the usual particle updates (insertions, dele¬ 
tions, displacements and rotations), combined with “bi¬ 
ased” insertions and deletions, which attempt to add (or 
remove) a colloid in a lock-key bound state, subject to 
satisfying detailed balance. This innovation, (described 
further in Appendix B), increases the efficiency of the 
simulation by up to 4 orders of magnitude. 

Fig. 2(a) summarises the phase behavior of the in¬ 
dented colloids. The state point of the system is specified 
by three parameters (p, £lK)£bb) where p is the number 
density of colloids. The solid line in Fig. 2(a) indicates 
values of (elKiEbb) for which a liquid-vapour critical 
point exists, at some critical density p“. For values of 
(elKjEbb) above this line, there exist values of the den¬ 
sity for which liquid-vapor phase coexistence occurs: see 


for example Figs. 2(b,c). The structure of these phases is 
discussed further below. Dotted lines in Fig. 2(a) show 
how the parameters of the effective model are related 
to depletant parameters {q,r]s)- for smaller q, lock-key 
binding is more favorable than back-to-back binding. In¬ 
creasing 77s at constant q leads to an increase in both elk 
and Ebb, along the dotted lines. To assess the accuracy 
of the effective potential, we used the GGA within a re¬ 
stricted Gibbs ensemble [31] to locate the critical point in 
a system of indented colloids with explicit depletant, for 
q = 0.2. We find that the critical parameters lie within 
10% of the results obtained with the effective potential. 
However, with explicit depletant, it was not possible to 
investigate the structure of liquid and vapor phases, nor 
the behavior for smaller q, due to the computational cost 
required. Hence our use of the effective potential in this 
work. 

The binodal curves associated with liquid-vapor phase 
coexistence are shown in Fig. 2(d), as a function of the 
depletant volume fraction rjg, for three values of q. As 
the depletant particles get smaller (decreasing q at fixed 
77s), the attractive forces between colloids get stronger, so 
phase separation occurs for smaller values of rjs- However, 
the most striking effect in Fig. 2(d) is the strong decrease 
(by more than a factor of 2) in the critical density of the 
colloids p“. The origin of this effect is the very strong 
lock-key binding that occurs when q is small. 

Figs. 2(b,c) illustrate the unusual structures of the co- 























3 


existing phases for q = 0.126, corresponding to strong 
lock-key binding. There are large void spaces within the 
liquid: we identify these by inserting spherical “ghost” 
particles of size a into the system, wherever this is possi¬ 
ble without overlapping an existing colloid (ghost parti¬ 
cles may overlap with each other). These particles, shown 
in Fig. 2(c), highlight the existence of the voids. At this 
state point, the depletant volume fraction is 5% greater 
than its critical value so the system is significantly 
outside the critical regime: these voids are not associ¬ 
ated with critical fluctuations, but are intrinsic to the 
liquid states that occur for small q. 

To reinforce the unusual nature of the liquid state at 
small q, we note that liquid configurations for larger q 
lack these large voids: they resemble “normal” colloidal 
liquids of spheres. The largest void in a typical liquid 
configuration at q = 0.2 and rjs = l.OSryg is shown in 
Fig. 2(e), and is much smaller than the corresponding 
void in Fig. 2(f). The ghost particles described above 
allow this effect to be analysed quantitatively: the typical 
volumes associated with the voids are shown in Appendix 
C. We have also analysed the state points at which the 
ghost particles percolate the system, as shown by the 
“void percolation” lines in Fig. 2(d). On the low-density 
side of these lines, the voids form a percolating network, 
and we expect that a test particle of a size comparable 
with a colloid can travel freely through the liquid: this 
may be a useful criterion for identifying “empty” [23- 
25] or “porous” liquids. At higher densities, the liquid 
resembles more familiar colloidal liquids. The key point 
to be noted from Fig. 2(d) is that decreasing q leads to an 
increasing range of density (p) and interaction strength 
{r]s) for which porous liquids exist. 

The unusual structure of these liquid and vapor states 
originates from a hierarchy of energy scales - the sys¬ 
tem supports strong lock-key bonds, and weaker back- 
to-back binding. To illustrate this effect, imagine that 
the colloids form linear chains, connected by strong lock- 
key bonds [17]. For large Elk, almost all of the avail¬ 
able lock-sites are participating in binding - the weaker 
back-to-back binding then provides both inter-chain and 
intra-chain interactions, which can cause the chains to 
aggregate or collapse, as happens in solutions of poly¬ 
mers [20-22, 32]. However, chains of colloids linked by 
depletion forces are in fact akin to “living polymers” [33], 
in that they are continuously breaking and re-forming in 
the equilibrium state. One can associate the liquid state 
with an aggregated state of many polymers, while the 
unusual vapor phase contains polymeric chains that tend 
to collapse into compact states [recall Fig. 2(b)]. 

In fact. Fig. 2(d) closely resembles the behaviour of 
polymers in solution, interacting via a relatively weak 
non-specihc attraction [20-22]. As the polymer length 
increases in that system, the critical point moves to 
lower polymer density and weaker non-specific interac¬ 
tions. These systems also exhibit strong “field-mixing” 
at the critical point [29, 30]: the natural order parameter 
for the liquid-vapor phase transition is not the density of 



FIG. 3. (a) Network motifs for lock-key bonds. Circles in¬ 
dicate particles and arrows indicate lock-key bonds, directed 
from the “lock” to the “key”. The most important local en¬ 
vironments are linear chains (one outward and one inward 
bond) and branching points (two inward and one outward 
bond). The converse branching case (one inward and two 
outward bonds) is not possible, due to the particle shape, 
(b) Clusters of particles can be characterised by networks of 
arrows. Once a connected cluster of particles includes a single 
loop, it can grow further only by forming inward bonds (some 
possible locations for additional bonds are shown with dotted 
arrows). It follows that each connected cluster can include at 
most one loop. 


the system, but rather a linear combination of the density 
and the energy. We show evidence in Appendix B that 
the liquid-vapor critical points of the indented colloids 
are in the Ising universality class (as expected), but with 
strong held-mixing, as in the polymer case [21]. This lat¬ 
ter effect rehects the unusual structures of the liquid and 
vapor states, especially the strong asymmetry between 
the structures of liquid and vapour states in Fig. 2(b,c), 
(that is, particle-hole asymmetry). 

This analogy with polymers leads to a prediction for 
the colloidal system: in the limit of long polymers, one 
expects the liquid-vapor critical point to occur at the 
0-temperature of the polymer [20-22], where the chain 
statistics are those of a simple random walk. This re¬ 
quires a weak non-specihc attraction, to overcome the 
excluded volume interactions that cause the polymer to 
swell. For the colloidal system, the long-polymer limit 
corresponds to £lk — > c», in which case one would ex¬ 
pect the liquid-vapor critical point to occur when the 
non-specihc (back-to-back) interaction ebb balances the 
excluded volume swelling effect. This implies that the 
line of critical points in Fig. 2(a) should tend to a non¬ 
zero value of Ebb, as Elk —>■ oo, consistent with our re¬ 
sults. 

However, while this analogy between chains of indented 
colloids and linear living polymers is appealing, it misses 
an important feature of the colloidal system. For the 


4 


particle shape considered here, these colloids readily form 
branched chains. In patchy-particle models, the presence 
of branching is sufficient to drive a liquid-vapour phase 
transition [23-25], but this is not the case for lock-and- 
key colloids. The reason is the directed nature of the lock- 
key bond. In Fig. 3, we illustrate different local structures 
in these systems, and the construction of a bond network. 
Lock-key bonds are indicated by arrows that point from 
the “lock” particle to the “key”. In this representation, 
the crucial feature is that a particle may have several 
inward bonds, but at most one outward bond. It follows 
that a cluster of particles connected by lock-key bonds 
may contain at most one loop, as illustrated in Fig. 3(b). 

This feature has important implications for phase 
behaviour - a condensed liquid phase can be stable 
only in the presence of percolating clusters that con¬ 
tain large numbers of internal loops. This effect can 
be illustrated using Wertheim’s theory of associating flu¬ 
ids [17, 23, 24, 34], which divides particles into “quasis¬ 
pecies” according to their local bonding, and predicts the 
number densities associated with these species. For in¬ 
dented colloids, we define pl as the number density of 
unoccupied “lock sites” (colloidal indentations), and px 
as the number density of free “key sites” (concave colloid 
surfaces) to which no particle is bound. The theory pro¬ 
vides “mass-action” equations which relate the numbers 
of bonds in the system to pn and px [17]: see Appendix 
D. 

Within Wertheim’s theory, the pressure P of a fluid of 
clusters satifies [24] 

pp = Ptree + ^Prei (1) 

where ptree is an estimate of the number density of 
bonded clusters in the system, and AP^ef is a contri¬ 
bution coming from a hard-particle reference system. 
The density ptree counts the number of clusters under 
the assumption that they are all “tree-like” (contain¬ 
ing no loops). For indented colloids without back-to- 
back interactions, we have ptree = Pl (in the absence 
of loops, every cluster has exactly one unoccupied lock 
site). However, liquid-vapor phase transitions require 
(dP/dp) < 0; it is easily shown within Wertheim’s the¬ 
ory that (dpi^/dp) > 0, and we expect on very general 
grounds that APj-ef is increasing in p. Hence, in the ab¬ 
sence of back-to-back binding, the theory predicts that 
no phase transition is possible: see Appendix D. 

However, if we include back-to-back binding within the 
theory, we find ptree = Pl—Pkk, where pxx is the number 
density of back-to-back bonds. Both pL and pxx are 
increasing functions of p - in this case the theory predicts 
that phase transitions can occur only when ptree < 0. 
Note however that if ptree is negative then it no longer 


represents a good estimate of the number of clusters in 
the system: due to the construction of ptree, this signifies 
that a proliferation of loops is necessary for condensation 
into a liquid phase. 

The theory used here treats back-to-back binding in 
a schematic way and is therefore not quantitatively ac¬ 
curate. Nevertheless, this analysis illustrates the impor¬ 
tance of cluster loops in liquid-vapour phase transitions. 
By contrast, the phase transitions that occur in patchy 
colloids are driven by branching and do not require “non¬ 
specific” interactions such as the back-to-back binding 
considered here - this is possible because there is no 
difference between the “inward” and “outward” bonds 
shown in Fig. 3, so there is no constraint on the num¬ 
ber of loops within a single cluster. In those cases, the 
Wertheim theory is quantitatively accurate [23, 24]. 

These results for indented colloids illustrate the subtle 
role of the bonding topology in hierarchical fluids. We 
see that loops are essential for condensation - they may 
arise either from “undirected” branching (as in patchy 
colloids) or from non-specific interactions (as in polymer- 
solvent criticality and indented colloids). Yet another 
point of reference is provided by dipolar fluids [35-37], 
where formation of ring (or loop) structures acts to sup¬ 
press condensation, due to the absence of non-specific 
interactions between closed rings, and the low probabil¬ 
ity of branching. Thus, while the formation of colloidal 
polymers may appear similar for several kinds of parti¬ 
cle, the energy scales associated with liquid-vapor phase 
separation can be quite different, as can the properties 
of the coexisting phases. 

In conclusion, indented colloids with depletion inter¬ 
actions support unusual liquid and vapor phases, which 
arise from a hierarchy of energy scales for lock-key 
and back-to-back binding. When the energy scales are 
well-separated, these systems also share some properties 
with empty liquids [23-25], dipolar fluids [35-37], and 
polymer-solvent systems [20-22, 32], but we emphasise 
that these systems have different mechanisms of conden¬ 
sation, due to the different topological properties of the 
bond network in each case. For small q, the liquid phases 
of indented colloids are “porous”, being characterised by 
large percolating voids, the typical size of which can be 
tuned by varying the properties of the depletant fluid. 
These states resemble “equilibrium gels”, which are not 
dynamically arrested [23, 25]. The fine experimental con¬ 
trol that can be achieved by manipulation of particle 
shape and depletant properties [14] makes indented col¬ 
loids ideal systems for further study in this direction. 
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Appendix A: PARAMETERISATION OF 

EFFECTIVE POTENTIAL 

1. Theoretical motivation 

As described in the main text, we use an effective po¬ 
tential that captures the important features of the de¬ 
pletion interaction between lock-and-key colloids. The 
effective interaction potential is v{ri 2 ,ni,n 2 ) where the 
vector between two particles is ri 2 and their orientations 
are ni. 2 . Our goal is to parameterise the function v so 
that our simplified model accurately represents the be¬ 
haviour of a mixture of hard indented colloids and ideal 
depletant particles. The optimal effective model would 
match two-body interactions exactly: 

exp[-/3?;(ri2,ni,n2)] = 5o(j’i2, Wi, ^ 12 ) (Al) 

where go is the two-particle distribution function for the 
colloids, normalised so that go = 1 aX large distances. 
The subscript ‘0’ indicates that this function is calculated 
in a dilute limit, where the colloid density is taken to zero 
(at fixed depletant volume fraction rjs)- 

The function go can be obtained computationally, but 
it depends on the distance ri 2 = |t’i 2 | and on three differ¬ 
ent angles. Parameterising this non-trivial function of a 
four variables is rather difficult - we have used numerical 
simulation to identify the most important features of this 
function, which we incorporate into a simplified effective 
potential. 

The radial distribution function in the dilute limit may 
be obtained as 

, , f dnidn2 _ , , , . 

9o{r) = J 9o{ri2,ni,n2) (A2) 

where both integrals run over the entire unit sphere. 
Fig. 4 shows an example of this function, in a situation 
where both lock-key and back-to-back binding are signfi- 
cant. The two separate peaks indicate the different bind¬ 
ing mechanisms, but this representation does not reveal 
the orientation dependence of the interaction. 

Fig. 5 shows the choice of co-ordinate system that we 
have used to capture this dependence. In lock-key bind¬ 
ing, one particle plays the role of the “lock” and the other 
the “key”. The definition cos = max(ni-ri 2 , —n 2 -T'i 2 ) 
takes the smaller of the angles between ri 2 and the direc¬ 
tors ^ 1 , 2 ; the particle associated with this smaller angle 
plays the part of the lock. We define cos0(^ = dc/cr so 
that 9^ corresponds to angle between the director and a 
vector pointing out through the “lip” of the lock mouth. 
Lock-and-key binding happens only for < 9^ while 
back-to-back binding happens for fi'R > ^R- 

Since back-to-back binding takes place between two 
convex surfaces, this part of the interaction depends 
weakly on angular co-ordinates, so the effective potential 
includes a square well of range tbr = ITcr in that case. 
The contribution to the effective potential v is therefore 
ubb = —£bb0(6*r — ^r)0(1-1o' — ri 2 ) where 0(a;) is the 



FIG. 4. Radial distribution function for two indented colloids, 
with depletant size ratio, q — 0.2, and gs = 0.30. (Recall dc = 
0.6 (t.) A thin dashed line shows the corresponding function 
in the absence of depletant. 




FIG. 5. Go-ordinate system for interactions among indented 
colloids. The angles between the interparticle vector ri 2 and 
the directors ni ,2 are labelled as shown: the labels 6 'r,i are 
chosen such that 6b. < di. The angle between the directors, 
is (f>. 


step function. For lock-and-key bonding, the angular de¬ 
pendence is less simple, but the dominant effect can be 
accounted for through the angle which measures the 
rotation of the “key” particle within the mouth of the 
lock. Hence the interaction potential for lock-key binding 
is ulk = -£lkW (0)0(6'b - 6'R)0(dc + O.lcr - ri 2 ) where 
W{(p) is the piecewise linear function of cos p shown in 
Fig. 1 of the main text. In fact, this simple function 
provides an almost quantitative description of the p- 
dependence of the effective potential [38]. 


2. Numerical parameterisation 


Having specified the form of the effective potential, we 
now describe how the parameters (elKi^bb) are chosen. 
The (effective) second virial coefficient for the indented 
colloids in the presence of the depletant is [39] 


B2 = 2tt 



drr^]! - go{r)] 


(A3) 
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We split B 2 into two pieces, associated with lock-key and 
back-to-back binding: 


^2 ^ = 27r /" drr^[l-go{r)] 

Jo 

poo 

S™=27r/ dvr'll!-go (r)] (A4) 

J a 


These two parameters indicate the strength of the two 
kinds of binding. A representative plot of go(^) is shown 
in Fig. 4. The peak at r « dc is associated with lock-key 
binding, and its area sets the strength of this interac¬ 
tion. The peak at r « ct is associated with back-to-back 
binding. 

From (A2), and B^^ can be written as integrals 
of 1 — 5 o- Motivated by (Al), our procedure is therefore 
to choose Slk, £bb so that the corresponding integrals of 
1 — exp[—/ 3 '!;(ri 2 , ni, ^ 2 )] are equal to B 2 ^,Bf^. These 
integrals are contributions to the second virial coefficient 
of the effective model, so we denote them by S 2 and 

rBB 

^2,eS- 

Since the potential is a square well and the hard- 
particle interactions have range at most cr, it is easily 
verified that 


= I'm - nLK(l - ) 

Bf^e = VBB{l-e^^-) (A5) 

where accounts for the excluded volume between two 
indented colloids, and ulKj'Cbb are phase space volumes 
associated with the two binding modes. 


3. Computational implementation 


To estimate goii^) for the system with depletant, we use 
the Geometrical Cluster Algorithm (GCA) to simulate 
two colloidal particles in a box, mixed with depletant, as 
in [17] (this procedure accounts for finite-size corrections 
to g{r)). Representative results for B^^, Bf^ are shown 
in Fig. 6 , for several size ratios q. As expected, these pa¬ 
rameters become large and negative as the nanoparticle 
volume fraction 77 s increases. 

The corresponding quantities B^^g, were calcu¬ 

lated numerically by the same method, but using the ef¬ 
fective potential instead of an explicit depletant. Results 
are shown in Fig. 7, including fits to (A5), from which 
we estimate the small volume associated with lock-and- 
key binding to be ulk = (3.3 x Comparison 

of Figs 6 and 7 allows straightforward mapping from the 
depletion parameters (g, rjs) to the parameters (elki Ebb) 
of the effective potential. 
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FIG. 6. The contributions to the second virial coefficient for a 
system of two locks, interacting with a depletant of penetrable 
spheres. 


Appendix B: SIMULATION METHOD 


1. Biased insertions in the grand canonical 
ensemble 


The lock-key interactions in our system are typically 
strong, implying a strongly negative chemical potential 
and consequently a very low acceptance rate for stan¬ 
dard particle transfers (insertions and deletions) in grand 
canonical Monte Carlo (GCMC) simulations. To counter 
this we introduce an additional biased transfer update in 
addition to the standard one. This operates by prefer¬ 
entially attempting to insert an indented particles in a 
bound configuration around a target colloid. 

To implement the “biased insertion”, we first pick a 
target colloid randomly from the existing ones in the box. 
The location of the centre of the colloid to be inserted is 
then chosen randomly from within a narrow spherical 
shell surrounding the target. The shell radii are chosen 
to represent the typical separation of bound colloids. For 
a given point in the shell there will be a range of solid 
angle for the particle director which avoids hard particle 
overlaps, and we choose the director of the new parti¬ 
cle uniformly from within this range. The biasing factor 
needed to maintain detailed balance is therefore a prod¬ 
uct of the shell volume and the solid angle that avoids 
overlap. This biasing factor can also be calculated for a 







FIG. 7. Contributions [see (A4)] to the second virial coef¬ 
ficient for a system of two indented colloids interacting via 
the effective potential. Dashed lines are fits to (A5) with 
Wni = 1.94 (t^, ulk = (3.3 x and hbb = 0.44(t^ 


bound particle which is to be removed from the system, 
based on its separation from its bound neighbour, and so 
enters symmetrically in the grand canonical acceptance 
probabilities. 


2. Determination of binodals, critical points and 
field mixing 

Our GCMC simulations sample the fluctuations in the 
colloid number density p = N/V and energy density 
u = E/V. Well within the two-phase region, the prob¬ 
ability distribution P{p) exhibits a pair of widely sepa¬ 
rated peaks centered on the coexistence densities. Ac¬ 
cordingly p serves as an order parameter for the transi¬ 
tion, and the coexistence densities (binodal) can be sim¬ 
ply read off from the peak positions [29]. Closer to the 
liquid-vapor critical point, the fluctuations in p and u 
assume a universal character which can be probed via 
the properties of the joint distribution P{p,u) [29,30]. 
In models that exhibit particle hole symmetry, the order 
parameter distribution at criticality matches the (inde¬ 
pendently known) universal fixed point form of the Ising 
magnetisation distribution P*{m). Effecting this match¬ 
ing provides a means of estimating the critical parame¬ 
ters. However, since our model lacks particle-hole sym¬ 
metry, field mixing is expected to occur. This implies 



FIG. 8. (a) Critical point order parameter distributions P{x) 
with X = a{p — pc))- The figure shows data for q = 0.2 and 
q = 0.126 and the distributions are scaled to unit norm and 
variance via the choice of the parameter a. Also shown for 
comparison is the fixed point Ising magnetisation distribution 
P*{m). (b) The corresponding normalised distribution of the 
field mixed order parameter P{y) with y = b(A4 — A4c) and 
b chosen to scale the distributions to unit variance. Box sizes 
were £ = 12ct for the q = 0.126 system and £ = 7.5a for the 
q = 0.2 system. 


that the scaling equivalent of m is the linear combina¬ 
tion M. = p — su [29,30]. (Strictly speaking, a non-linear 
combination may be appropriate 40], but the linear com¬ 
bination is a simple and practical choice.). To estimate 
the critical parameters the field mixing parameter s must 
be tuned along with p, and T such as to best match P{A4) 
to P*{m) [29]. 

We have used these techniques to determine the critical 
points and coexistence binodals for a selection of values 
of q in the interval q = [0.25,0.1]; the results are shown 
in Fig. 2(a,d) of the main text. Our universal distribu¬ 
tion matching method also delivers estimates for the field 
mixing properties of the critical fluid. For q = 0.2, a sat¬ 
isfactory matching of the number density distribution to 
the Ising fixed point form can be obtained as shown in 
Fig. 8(a), so field mixing is negligible in this case. By 
contrast, for q — 0.126, no satisfactory matching can be 
obtained under the assumption s = 0. The situation is 
greatly improved, however, by forming the distribution 
of the mixed field ordering operator P{A4), as shown in 
Fig. 8(b), which corresponds to s = —0.058(1) 
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q = 0.126 


FIG. 9. (a) Fraction of space available to spherical 

“ghost particles” of size in liquid states at different 

size ratios q. (The probe diameter is measured in units of 
the colloid diameter a.) For the porous liquids at 

q = 0.126 have many more cavities where ghost particles will 
fit, compared with the dense liquids at g = 0.2. (b) Depiction 
of void spaces in representative configurations, for = 1. 

Each connected void space is shown with a different color, to 
illustrate their sizes and connectivity. 


Appendix C: VOID SIZES IN LIQUID STATES 

In Fig. 9, we show the fraction of the total vol¬ 
ume of the system that is accessible to “probe particles” 
(or “ghost” particles) of size These calculations 

were performed at states where ijs = l.OSr^g. The porous 
structure of the liquid at g = 0.126 means that is 
greater by several orders of magnitude than the corre¬ 
sponding state a.t q — 0.2. We also show the void spaces 
associated with typical configuration, at = 1. Note 

that the sizes of the typical voids increase along with 
their volume fraction; at q = 0.126, the largest cluster 
has percolated the system. 


Appendix D: WERTHEIM THEORY 


tion, where the density is decomposed as a sum of partial 
densities that describe particles in different local environ¬ 
ments. The second part is an (approximate) thermody¬ 
namic perturbation theory (TPT) which facilitates ex¬ 
plicit calculations. In [17], we developed this theory for 
a system of indented colloids that forms linear chains, 
but both chain branching and back-to-back binding were 
neglected. 

The generalisation of this theory to include these ef¬ 
fects is straightforward, although some book-keeping is 
required to keep track of the various local environments. 
Briefly, the interaction between particles is broken down 
into three parts, which take care of particle overlaps, lock- 
and-key binding, and back-to-back interactions. This de¬ 
composition is performed at the level of the Mayer /- 
function, so the diagrammatic expansion of the free en¬ 
ergy contains three kinds of interaction. 

Within this diagrammatic expansion, Wertheim’s the¬ 
ory requires that particles be classified according to their 
local bonding. In practice, we keep track of the number 
of inward and and outward lock-and-key bonds, and the 
number of back-to-back bonds. The constraint (recall 
Fig. 3 of main text) that each particle has at most one 
outward bond is enforced at the level of the diagrammatic 
series (this is the key advantage of Wertheim’s approach). 


1. Quasispecies densities and mass-action equations 

To develop the theory, let pijk be the number of par¬ 
ticle with i inward bonds, j outward bonds and k back- 
to-back bonds. These are the quasispecies densities of 
Wertheim’s theory [34]. For simplicity we assume that 
particles may have at most two inward lock-and-key 
bonds {k < 2), at most one outward lock-and-key bond 
(j < 1) and at most one back-to-back bond {k < 1). 
Further, we assume that that particles with back-to-back 
bonds have no inward lock-and-key bonds and vice versa, 
so that only one of j or k may be greater than zero. Even 
with these constraints, there are still eight distinct local 
environments, but the advantage of the TPT is that the 
various pijk have fairly simple relations between them. 

We postpone the detailed analysis to a later work and 
quote only the main results, with their physical inter¬ 
pretation. We use the notation po = pooo, for brevity. 
“Chemical equilibrium” between free monomers and par¬ 
ticles with a single inward bond (and no other bonds) 
means that 


Poio = Kpopi^ (Dl) 


where 


In the main text, arguments based on Wertheim’s the¬ 
ory of associating fluids [34] were used to investigate the 
effects of branching and non-specific binding on phase 
behaviour. The theory consists of two main parts. The 
first is an exact representation of a fluid’s partition func- 


Pl — Po + Pwo + P 020 + Pool (D2) 

is the density of unoccupied lock sites (see also the main 
text), and K is an equilbrium constant which, at the sim¬ 
plest level, is equal to More generally, K may be 
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obtained within the TPT as an integral involving the ra¬ 
dial distribution function g{r) in a hard-particle reference 
system. 

We also define 

Pk = Po + Poio 
Plk = Poio + Piio 

PhhK = P020 + Pl20 (D3) 

which correspond to number densities of key surfaces 
with 0, 1, and 2 inward bonds. Also, 

Pkk = 5 (pool + Pioi) (D4) 

is the number density of back-to-back bonds. 

Chemical equilibria between other particle quasispecies 
yield 

Pool = ATbbPoPk 

P020 = oAT^poPl (D5) 

where Kbb is an equilibrium constant for back-to-back 
binding, and a is a geometrical factor that accounts for 
the reduced volume available when binding a second lock 
onto a key where a lock is already bound. (The theory 
provides integral expressions for these quantities; in gen¬ 
eral we expect a to depend weakly on the energy scales 
Elk, £bb while Kbb oc The Wertheim theory also 

builds in independence of bonds for the lock- and key- 
sites of a given particle. This means for example that 
Piio = PiooPoio/po, and similarly pi 2 o = PiooPo 2 o/po and 
Pioi = PiooPooi/po- 

Combining all the definitions and theoretical results 
given so far, we arrive at 

Plk = KpbPk 
Pllk = clKpbPk 

Pkk = A'bbPk (D6) 

which are the “mass-action” equations described in the 
main text. 

Since the total density p is equal to the sum of all the 
partial densities po,Pooi etc, we also have 

P = Pk + Plk + Pllk + 2pkk (D7) 

which simply means that all key-sites must be either free 
or involved in some kind of bonding. Combining (D1,D5) 
yields piooPL = (poio + 2po2o)PK from which one may 
derive 

P = Pl + Plk + 2pllk (D8) 

This equation enforces the constraint that all lock-sites 
are either occupied or unoccupied, analogous to (D7). 


2. Pressure 

Following [23], we obtain the pressure of the system 
within the TPT as a sum of two terms 

pp = Ptree + ^Prei (D9) 

where ptree is an estimate of the number of clusters in 
the system and AP^ef = PPrei — p is a contribution from 
the hard-particle reference system, which has a reduced 
pressure pPref that depends on the total density p but 
not on the quasispecies densities. As discussed in the 
main text, ptree is equal to the number of clusters in 
the system if there are no internal loops within clusters 
(in that case, formation of any bond decreases the total 
number of clusters by 1). This “tree-like” assumption 
motivates the notation ptree- 

The theory described here gives 

Ptree = PL — PKK- (DIO) 

To see this, note that if all clusters are tree-like then 
their number is equal to the difference between the total 
number of particles and the total number of bonds (either 
lock-key or back-to-back). The number of lock-key bonds 
is Plk + 2plkk = P — Pl and the number of back-to-back 
bonds is pkk, yielding (DIO). 

a. Conditions for phase separation 

Liquid vapour phase transitions occur if 
{dP/dV)N,T > 0 for some volume V. Here the 
derivative is taken at a constant value of the total 
particle number JV, but the numbers of particles 
within each quasispecies depend on V. The reference 
contribution APref is decreasing in D so a necessary 
condition for phase separation within Wertheim’s 
theory is that {dptree/dV)N,T > 0. At constant 
T, Ptree is a simple function of p, which implies 
y(pPtreeldV)M,T = ~p(dptree/dp), SO phase Separation 
can occur only if 

(dp 

tree /dp) < 0. (Dll) 

(More precisely, one requires dptree/dp < 

-d(APref)/dp.) 

The relation (Dll) provides a useful necessary condi¬ 
tion for phase separation. In particular, for the system 
considered here, the densities pl, Pk, Pkk, etc are all 
monotonic in p (see Sec. D 2 b below). In the absence of 
back-to-back binding, one has ptree = Pl from (DIO), the 
pressure is monotonic in p, so the theory predicts that 
no phase transition is possible. 

Physically, the idea is that internal loops within 
bonded clusters of particles are suppressed by the di¬ 
rected lock-key bond and the fact that a lock-site can 
bind to at most one key-site. The Wertheim theory cap¬ 
tures this property. 
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b. Monotonicity of pi,, pK 

To complete the argument that phase transitions 
within the theory requires it'sH > 0, we prove that pn 
and pk depend monotonically on p. Starting from (D8), 
we use (D6) to express all quantities in terms of Pl,Pk 7 
and we consider small variations in (p, pl, Pk), obtaining 


and again considering small variations in pu, pk yields 
(5pL(l+2aitr^PLPK) = <5 pk( 1—ait'^PL+ditiBBPK)- (D14) 
Using (D13) again gives 

dpL(l + 2aif2pLPK) = —(^PL + 2 KbbPk). (D15) 

PK 


Sp = dpL(l + KpK + iaK'^phPK) + SpKiKpL + 2aK‘^p^) 

(D12) 


From (D7,D8), we have 


Pl + Pllk = Pk + 2pKK (D13) 


It follows from (D12,D15) that Sp and SpK are both 
monotonic in dpn, which implies that (pl,Pk) are both 
monotonic in p, as are (pkk, Plk, Plkk)- Given that 
Ptree = PL — PkKj it follows that ptree must be mono- 
tonic if Pkk = 0, in which case (dptree/dp) > 0 and the 
theory predicts no phase transition. 



